Technical issues in mixed models

During this computer lab we will 1) Try to reproduce some results from the second lecture. 2) Revisit examples from earlier in the mixed models part, and do some proper model building and checking. 3) Expand on the Reisby analyses using polynomial functions for time. 4) Analyze data from a multi-center, longitudinal trial. 5) (Optional) Examine effects of centering explanatory variables. 6) (Optional) Take a look at a (two-level) dataset with complicated residuals.

Note: Exercise 3 is intended to help you familiarize yourself with mixed models in R. If you are already comfortable with the software, feel free to skip it.

Before we get started, we’ll first load a few packages we will need for our analysis.

library(nlme)
library(car)
library(ggplot2)
library(lmerTest)

Exercise 1

Repeat the analysis of the Reisby data in R.

load(file = "reisbywide.Rdata")
load(file = "reisbylong.Rdata")
reisby.long$id <- factor(reisby.long$id)
reisby.long$time <- reisby.long$week
reisby.long$week <- factor(reisby.long$week)
  1. Modelling time (linear), treatment and time*treatment in the fixed part of the model, decide on the best random structure for the data (random intercepts and/or slopes).

We’ll start by modelling time as continuous (and linear) and examine a few linear mixed effects models. All models include fixed effects for time, endo, and their interaction. The first includes only a random intercept per patient, and the second a random intercept & random slope for time.

# mixed model with random intercept, fixed effects for time(linear!), endo & interaction
lme.ril<-lme(fixed=hdrs ~ time*endo, random=~1|id, data=reisby.long, na.action="na.omit", method="ML")
# mixed model with random intercept & random slope for time(linear), fixed effects for time, endo & interaction
lme.ris<-update(lme.ril, random=~time|id)
anova(lme.ris, lme.ril)
##         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## lme.ris     1  8 2230.929 2262.345 -1107.465                        
## lme.ril     2  6 2294.137 2317.699 -1141.069 1 vs 2 67.20798  <.0001

Not surprisingly, the model with a random intercept and random slope is significantly better than the model with just a random intercept per person.

  1. Using the best random part from part (a), see if you can reduce the fixed effects in the model.
lme2.ris<-update(lme.ris, fixed=hdrs ~ time+endo)
lme3.ris<-update(lme.ris, fixed=hdrs ~ time)
anova(lme.ris, lme2.ris, lme3.ris)
##          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## lme.ris      1  8 2230.929 2262.345 -1107.465                        
## lme2.ris     2  7 2228.933 2256.422 -1107.467 1 vs 2 0.004160  0.9486
## lme3.ris     3  6 2231.037 2254.599 -1109.519 2 vs 3 4.104108  0.0428

We’ll re-fit the final model now using restricted maximum likelihood estimation to get proper variance estimates. And we’ll ask for approximate 95% confidences intervals for the various parameters in the model:

# Best model is lme3.ris.CAR1; re-run it now with REML
lme3.ris.reml <- update(lme3.ris, method="REML")
summary(lme3.ris.reml)
## Linear mixed-effects model fit by REML
##  Data: reisby.long 
##        AIC      BIC    logLik
##   2231.918 2255.447 -1109.959
## 
## Random effects:
##  Formula: ~time | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 3.597855 (Intr)
## time        1.458071 -0.282
## Residual    3.494655       
## 
## Fixed effects: hdrs ~ time 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 23.577044 0.5499025 308  42.87495       0
## time        -2.377047 0.2103445 308 -11.30073       0
##  Correlation: 
##      (Intr)
## time -0.45 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.73577812 -0.49824400  0.03320466  0.51450770  3.67933067 
## 
## Number of Observations: 375
## Number of Groups: 66
intervals(lme3.ris.reml)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                 lower      est.     upper
## (Intercept) 22.495003 23.577044 24.659085
## time        -2.790941 -2.377047 -1.963153
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: id 
##                            lower       est.      upper
## sd((Intercept))        2.7386421  3.5978545 4.72663347
## sd(time)               1.1428595  1.4580709 1.86022054
## cor((Intercept),time) -0.5608171 -0.2815758 0.05518275
## 
##  Within-group standard error:
##    lower     est.    upper 
## 3.194711 3.494655 3.822760

We found no evidence in our sample for different patterns over time for people with endogenous vs. exogenous depression (interaction not statistically significant), and no evidence in our sample that depression scores differ between people with endogenous vs. exogenous depression at any given time point. However, we do have evidence for decreasing depression scores over the six weeks. Patients treated with imipramine have an average weekly decline in HDRS of 2.4 points, 95% CI [2.0; 2.8]. There is considerable variation in the HDRS scores at baseline (standard deviation of the random intercepts is 3.6), and also considerable variation in the slopes of HDRS over time (standard deviation of the random slopes 1.5).

  1. (optional) Try fitting a few covariance pattern models, using time as categorical in the fixed part of the model.
# CPM, correlation = AR(1) homogeneous
gls.AR1 <- gls(hdrs ~ week*endo, correlation=corAR1(form = ~ 1 | id), data=reisby.long, na.action="na.omit", method="ML")
# corMatrix(gls.AR1$modelStruct$corStruct)[[1]]

# CPM, correlation = AR(1) heterogeneous
gls.hAR1 <- gls(hdrs ~ week*endo, correlation=corAR1(form = ~ 1 | id), weight = varIdent(form = ~ 1 | week), data=reisby.long, na.action="na.omit", method="ML")
# corMatrix(gls.hAR1$modelStruct$corStruct)[[1]]

# CPM, correlation = COMPLETELY unstructured
gls.unstruct <- gls(hdrs ~ week*endo, correlation=corSymm(form = ~ 1 | id), weight = varIdent(form = ~ 1 | week), data=reisby.long, na.action="na.omit", method="ML")
# summary(gls.unstruct)
# corMatrix(gls.unstruct$modelStruct$corStruct)[[1]]

To get the -2LL and the AIC of the homogeneous AR(1) model:

-2*logLik(gls.AR1)[1]
## [1] 2221.847
AIC(gls.AR1)
## [1] 2249.847

Now we’ll get the -2LL and AIC of the 3 different correlation structure models and put them into a nice table for comparison purposes:

fittable <- matrix(nrow=3, ncol=2)
fittable[, 1] <- c(-2*logLik(gls.unstruct)[1], -2*logLik(gls.AR1)[1],
                   -2*logLik(gls.hAR1)[1]  )
fittable[, 2] <- c(AIC(gls.unstruct), AIC(gls.AR1), AIC(gls.hAR1))
colnames(fittable) <- c("-2*logLik", "AIC")
rownames(fittable) <- c("gls.unstr", "gls.AR1", "gls.hAR1")
fittable
##           -2*logLik      AIC
## gls.unstr  2183.227 2249.227
## gls.AR1    2221.847 2249.847
## gls.hAR1   2207.462 2245.462

From the above table, we see that the AIC lowest is for the heterogeneous AR(1) model. We can also use the -2LL to compare nested models with one another using the LRT:

anova(gls.AR1, gls.hAR1)
##          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## gls.AR1      1 14 2249.847 2304.824 -1110.924                        
## gls.hAR1     2 19 2245.462 2320.074 -1103.731 1 vs 2 14.38499  0.0133

The LRT test indicates that heterogeneous AR(1) is statistically significantly better than homogeneous AR(1), p = 0.0133. The AIC would also agree in this case. Since the three models are not all nested, it’s better to use AIC for all model comparisons.

  1. Check the model assumptions of the final model, and interpret the results.

First, we’ll look at a Q-Q plot of the level 1 (child-level) residuals. Note that the code here is quite simple: in the nlme package, the function qqnorm() has been programmed to “communicate” with lme objects, and we only have to ask for a Q-Q plot of the object (model) in order to get a Q-Q plot of the standardized residuals.

qqnorm(lme3.ris.reml)

To check possible heterscedasticity and/or problems with linearity, we look at a scatterplot of the standardized residuals against the fitted values. Again, the plot() function works with lme objects, so we only need to plot the object to get the graph we want.

plot(lme3.ris.reml)

We can also examine the standardized residuals versus fitted values separately by endo/exo:

plot(lme3.ris.reml, resid(., type = "p") ~ fitted(.) | endo, abline = 0)

What do you see in these graphs? To get a sense of possible outliers, we can also look at boxplots of the residuals per subject:

plot(lme3.ris.reml, id ~ resid(.))

We see that some individuals have relatively little spread in their residuals, but most have residuals ranging from -5 or more to +5 or more. We also see that, for instance, subject 607 (4th subect from the top) has one very low residual, so one HDRS measurement that does not fit his/her predicted trajectory over time very well. The next subject down, 606, has one very high residual. You could go back to the original data and examine these cases in more detail if were concerned about particular outlier.

We can check poor fit of the model per school by plotting the observed versus fitted values per subject:

plot(lme3.ris.reml, hdrs ~ fitted(.) | id, abline = c(0,1))

The lines indicate perfect predictions, i.e. perfect agreement between observed (Y axis) and fitted (X axis). For most subjects, the values estimated from the model fit the original HDRS scores fairly well. However, we do see that some subjects have quite a bit of variation around the line of “perfect fit”. And that a couple of subjects (322, 328) have very little agreement between their observed and fitted HDRS.

Finally, we’ll get dotplots of the random effects (intercepts & slopes for time) per subject, to check for outliers. Again, the plot() function has been modified; when we apply it to the random effects of our model, we automatically get dotplots per subject. Note that we are not looking for a trend in this plot, simply for subjects with very low or high estimated random intercepts or slopes.

plot(ranef(lme3.ris.reml))

As we saw in the slides, there are a few subjects with higher than expected random slopes for time. Otherwise, there don’t appear to be any obvious outliers.

Exercise 2

We will re-analyze the schools data (school.Rdata) in R according to the strategy described in the lecture: starting with the full fixed model we first test the random effects, then the fixed effects.

load(file = "school.Rdata")
  1. Making use of the full fixed model from Part 1, test (with the likelihood ratio test) whether a random slope is necessary, or whether a random intercept model would be sufficient. Which model do you prefer?
### Model comparisons
# mixed model with random intercept & random slope, plus gender, school gender & school avg
sch.lme.1 <- lme(normexam~standlrt + factor(gender)+ factor(schgend) + factor(schav),
                 random=~standlrt | school, data=london, method="ML")

# mixed model with random intercept, plus gender, school gender & school avg
sch.lme.2 <- update(sch.lme.1, random=~1 | school)
# test random slope in full fixed model
anova(sch.lme.1, sch.lme.2)
##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## sch.lme.1     1 11 9300.414 9369.809 -4639.207                        
## sch.lme.2     2  9 9336.279 9393.057 -4659.139 1 vs 2 39.86525  <.0001
# random slope fits much better!

The LRT is highly significant (or: the AIC is considerably lower for the model with random slopes for standlrt), indicating that the model with a random slope fits much better. We cannot reduce the random effects in the model, so we stop here.

  1. Using maximum likelihood estimation and the “random part” chosen in part (a), try to reduce the fixed part of the model by removing non-significant explanatory variables.
# reduce fixed effects?
drop1(sch.lme.1, test="Chisq")
## Single term deletions
## 
## Model:
## normexam ~ standlrt + factor(gender) + factor(schgend) + factor(schav)
##                 Df    AIC     LRT  Pr(>Chi)    
## <none>             9300.4                      
## standlrt         1 9454.8 156.383 < 2.2e-16 ***
## factor(gender)   1 9322.7  24.320 8.159e-07 ***
## factor(schgend)  2 9302.3   5.927   0.05163 .  
## factor(schav)    2 9299.1   2.677   0.26227    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The variable schavg can be dropped. We’ll re-fit the model without it, and see if anything else can be dropped:

# schavg can be dropped
sch.lme.3 <- update(sch.lme.1, fixed= normexam~standlrt + factor(gender)+ factor(schgend))
drop1(sch.lme.3, test="Chisq")
## Single term deletions
## 
## Model:
## normexam ~ standlrt + factor(gender) + factor(schgend)
##                 Df    AIC     LRT  Pr(>Chi)    
## <none>             9299.1                      
## standlrt         1 9455.4 158.312 < 2.2e-16 ***
## factor(gender)   1 9321.7  24.659 6.844e-07 ***
## factor(schgend)  2 9301.4   6.267   0.04356 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# all remaining vars are significant

All remaining variables contribute significantly to better model fit, so we do no more further model reduction. Our final model is then the one with standardized London Reading Test, child’s gender and the school gender (all boys, all girls, or mixed).

  1. Once you have a final model, run it one last time using REML estimation for correct parameter estimates. Interpret the results of this model.
# Re-run final model with REML estimation
sch.lme.fin <- lme(normexam~standlrt + factor(gender)+ factor(schgend),
                   random=~standlrt | school, data=london, method="REML")
summary(sch.lme.fin)
## Linear mixed-effects model fit by REML
##  Data: london 
##        AIC      BIC    logLik
##   9321.058 9377.825 -4651.529
## 
## Random effects:
##  Formula: ~standlrt | school
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.2892361 (Intr)
## standlrt    0.1227839 0.574 
## Residual    0.7418142       
## 
## Fixed effects: normexam ~ standlrt + factor(gender) + factor(schgend) 
##                       Value  Std.Error   DF   t-value p-value
## (Intercept)      -0.1888921 0.05248756 3992 -3.598797  0.0003
## standlrt          0.5539994 0.02011862 3992 27.536647  0.0000
## factor(gender)1   0.1685081 0.03383955 3992  4.979619  0.0000
## factor(schgend)2  0.1797390 0.10170944   62  1.767181  0.0821
## factor(schgend)3  0.1744336 0.08078764   62  2.159162  0.0347
##  Correlation: 
##                  (Intr) stndlr fct()1 fct()2
## standlrt          0.313                     
## factor(gender)1  -0.302 -0.038              
## factor(schgend)2 -0.464  0.006  0.150       
## factor(schgend)3 -0.458  0.019 -0.230  0.240
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.84334881 -0.63773954  0.02664107  0.67776175  3.42763339 
## 
## Number of Observations: 4059
## Number of Groups: 65
intervals(sch.lme.fin)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                        lower       est.       upper
## (Intercept)      -0.29179702 -0.1888921 -0.08598716
## standlrt          0.51455566  0.5539994  0.59344313
## factor(gender)1   0.10216367  0.1685081  0.23485251
## factor(schgend)2 -0.02357532  0.1797390  0.38305323
## factor(schgend)3  0.01294139  0.1744336  0.33592577
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: school 
##                                lower      est.     upper
## sd((Intercept))           0.23580647 0.2892361 0.3547719
## sd(standlrt)              0.09046099 0.1227839 0.1666562
## cor((Intercept),standlrt) 0.23304066 0.5736876 0.7889340
## 
##  Within-group standard error:
##     lower      est.     upper 
## 0.7255837 0.7418142 0.7584078

Standardized LRT, child’s gender, and school gender all have a significant effect on the exam scores. The intercept is for boys (gender = 0) at mixed schools (schgend=1) with an average LRT (standlrt = 0). On average is the normalized exam score 0.55 SD higher when the standardized LRT is increased by 1 SD. Girls score, on average, 0.17 SD higher than boys, and children at single-sex schools score 0.18 (boys’ schools) or 0.17 (girls’ schools) SD higher than children at mixed-gender schools.

  1. Check the model assumptions for the final model: are the residuals, random intercept (and random slope if you used it) (approximately) normally distributed? (Note that some of the graphs for checking model assumptions are not easily produced in SPSS; see the extra slides on Moodle for help.)
### Model diagnostics
# Plots of residuals
qqnorm(sch.lme.fin)

plot(sch.lme.fin)

# Dotplot of residuals by school
plot(sch.lme.fin, school ~ resid(.))

# observed versus fitted values by school
plot(sch.lme.fin, normexam ~ fitted(.) | school, abline = c(0,1))

# Plot of random intercepts & slopes
plot(ranef(sch.lme.fin))

which.max(ranef(sch.lme.fin)[,2])
## [1] 53
ranef(sch.lme.fin)[53,]
##    (Intercept)  standlrt
## 53   0.5040069 0.3349675

No obvious problems with normality of residuals, distribution of residuals across fitted values. No one school stands out with high/low residuals. One clearly higher slope for stdlrt (school #53).

Exercise 3

In this exercise we will consider a possible quadratic effect of time in the Reisby example.

  1. Based on the spaghetti plots of the data (see the spaghetti plots in assignments part 2 exercise 1) one might suspect it is too simplistic to assume that the change across time is linear. Perhaps the trend in time is curvilinear? Test (with the likelihood ratio test) whether adding a fixed quadratic time parameter would improve the model. (Note: when using a quadratic term for time along with the linear term, it is better to first center time.)
# Center time variable and examine effect on random intercepts and slopes
# For simplicity: use model with fixed time effect, random effects for intercept + slope
reisby.long$ctime = reisby.long$time - 2.5
timenoc.mod1 <- lme(fixed=hdrs ~ time, random=~time|id, data=reisby.long, na.action="na.omit", method="REML")
timec.mod1 <- lme(fixed=hdrs ~ ctime, random=~ctime|id, data=reisby.long, na.action="na.omit", method="REML")
summary(timenoc.mod1)
## Linear mixed-effects model fit by REML
##  Data: reisby.long 
##        AIC      BIC    logLik
##   2231.918 2255.447 -1109.959
## 
## Random effects:
##  Formula: ~time | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 3.597855 (Intr)
## time        1.458071 -0.282
## Residual    3.494655       
## 
## Fixed effects: hdrs ~ time 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 23.577044 0.5499025 308  42.87495       0
## time        -2.377047 0.2103445 308 -11.30073       0
##  Correlation: 
##      (Intr)
## time -0.45 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.73577812 -0.49824400  0.03320466  0.51450770  3.67933067 
## 
## Number of Observations: 375
## Number of Groups: 66
summary(timec.mod1)
## Linear mixed-effects model fit by REML
##  Data: reisby.long 
##        AIC      BIC    logLik
##   2231.918 2255.447 -1109.959
## 
## Random effects:
##  Formula: ~ctime | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev   Corr  
## (Intercept) 4.341249 (Intr)
## ctime       1.458085 0.606 
## Residual    3.494638       
## 
## Fixed effects: hdrs ~ ctime 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 17.634426 0.5647201 308  31.22684       0
## ctime       -2.377047 0.2103457 308 -11.30067       0
##  Correlation: 
##       (Intr)
## ctime 0.493 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.73576223 -0.49824034  0.03320365  0.51449306  3.67933700 
## 
## Number of Observations: 375
## Number of Groups: 66

Observe that centering does not change estimates, Std. Error, and t-vales for the regression coefficient of time.

# Examine quadratic terms
# To avoid multicollinearity, it's actually better to use centered version of time when squaring!
# (Better still: orthogonal polynomials, not covered here)
# For simplicity: add ctime & ctime squared to model with fixed time effect, random effects for intercept + slope
reisby.long$ctimesq <- reisby.long$ctime^2
mod1 <- lme(fixed=hdrs ~ ctime, random=~ctime|id, data=reisby.long, na.action="na.omit", method="ML")
mod2 <- lme(fixed=hdrs ~ ctime + ctimesq, random=~ctime|id, data=reisby.long, na.action="na.omit", method="ML")
anova(mod1, mod2)
##      Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## mod1     1  6 2231.037 2254.599 -1109.519                         
## mod2     2  7 2232.626 2260.115 -1109.313 1 vs 2 0.4114832  0.5212

On the basis of either the AIC or LRT, we would conclude that model 2 (adding fixed quadratic effect for time) is no better than model 1 (with just linear effect).

  1. Would the model be improved if we also postulate that the quadratic time parameter individually deviates from average (i.e. is also random)?
mod3 <- lme(fixed=hdrs ~ ctime + ctimesq, random=~ctime + ctimesq|id, data=reisby.long, na.action="na.omit", method="ML")
anova(mod1, mod2, mod3)
##      Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## mod1     1  6 2231.037 2254.599 -1109.519                         
## mod2     2  7 2232.626 2260.115 -1109.313 1 vs 2  0.411483  0.5212
## mod3     3 10 2227.648 2266.917 -1103.824 2 vs 3 10.978143  0.0118

On the basis of either the AIC or LRT, we would conclude that model 2 (adding fixed quadratic effect for time) is no better than model 1 (with just linear effect). However, when we add a quadratic random term, the model fit is improved.

  1. The results of this last model suggest that although the trend across time is essentially linear at the population level, it is curvilinear at the individual level. Can you explain this phenomenon?

There are probably equal numbers of people with a convex and concave curve. On average, we have a linear time trend, though some individuals may have quadratic time trends.

Exercise 4

Last week we looked at a multi-center, randomized, double-blind clinical trial to compare three treatments for hypertension. One treatment was a new drug (A = Carvedilol) and the other two were standard drugs for controlling hypertension (B = Nifedipine, C = Atenolol). Twenty-nine centers participated in the trial and patients were randomized in order of entry. One pre-randomization and four post-treatment visits (at weeks 3, 5, 7 and 9) were made. We would like to see if there is a difference among the three treatments. The data can be found in the file bp.Rdata. Read the data into R. The research question is which of the medicines (treat) is more effective in reducing DBP. Since baseline (pre-randomization) DBP (dbp) will likely be associated with post-treatment DBP, we wish to include it here as a covariate.

Detach nlme library, attach lme4 (easier for fitting 3-level models) (You can leave them both attached, but there is some overlap in functions so it’s better to only have one open at a time)

detach("package:nlme")
library(lme4)
load(file = "bp.Rdata")
  1. First fit a two-level model, examining the effects of treatment, time and their interaction, adjusting for baseline dbp. Adjust only for multiple measurements per person (i.e. ignore clustering in centers) by including a random intercept and random slope for time per patient. Use ML estimation.
# NOTE: below I use the lmer function from the package lme4, because I think it is easier
# to define the random effects for the three-level model
# Two-level model (ignoring center): fixed effects time, tx, baseline DBP, time*tx
# random effects: intercept & time per patient
bp.2l.lme <- lmer(dbp ~ visit + factor(treat)+ visit*factor(treat) + dbp1 + (visit|patient), data=bp, REML=F)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00325065 (tol = 0.002, component 1)
summary(bp.2l.lme)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: dbp ~ visit + factor(treat) + visit * factor(treat) + dbp1 +  
##     (visit | patient)
##    Data: bp
## 
##      AIC      BIC   logLik deviance df.resid 
##   7501.6   7556.6  -3739.8   7479.6     1081 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9424 -0.5186 -0.0414  0.5045  3.4140 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  patient  (Intercept) 50.9908  7.1408        
##           visit        0.4885  0.6989   -0.47
##  Residual             33.0284  5.7470        
## Number of obs: 1092, groups:  patient, 288
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            44.82497    9.02458 294.45361   4.967 1.15e-06 ***
## visit                  -0.58399    0.15629 275.96392  -3.737 0.000227 ***
## factor(treat)B         -1.57395    1.59435 286.13525  -0.987 0.324377    
## factor(treat)C         -3.28604    1.58231 283.71256  -2.077 0.038726 *  
## dbp1                    0.50026    0.08692 289.38929   5.755 2.20e-08 ***
## visit:factor(treat)B    0.05733    0.22140 272.27807   0.259 0.795877    
## visit:factor(treat)C    0.04369    0.21842 270.12596   0.200 0.841599    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) visit  fct()B fct()C dbp1   vs:()B
## visit       -0.094                                   
## factr(trt)B -0.120  0.538                            
## factr(trt)C -0.078  0.542  0.490                     
## dbp1        -0.992 -0.001  0.035 -0.008              
## vst:fctr()B  0.059 -0.706 -0.770 -0.382  0.008       
## vst:fctr()C  0.063 -0.716 -0.385 -0.771  0.005  0.505
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00325065 (tol = 0.002, component 1)
  1. Now fit a three-level model, examining the fixed effects of treatment, time, their interaction, and baseline dbp, including a random intercept and random slope per patient and a random intercept per center. Use ML estimation. Note: it may be easier to switch the the lmer function in the lme4 package, allowing you to add random effects per level using + (effect|level) to the equation. If you have problems with convergence of a model, try adding the options control = lmerControl(optimizer = “optimx”, calc.derivs = FALSE, optCtrl = list(method = “nlminb”, starttests = FALSE, kkt = FALSE)) to your lmer model. For this to work you need to have the package optimx installed!
# Three-level model: fixed effects time, tx, baseline DBP, time*tx
# random effects: intercept & time per patient AND intercept per center
bp.3l.lme <- lmer(dbp ~ visit + factor(treat)+ visit*factor(treat) + dbp1 +
                    (visit|patient) + (1|center), data=bp, REML=F)
summary(bp.3l.lme)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: dbp ~ visit + factor(treat) + visit * factor(treat) + dbp1 +  
##     (visit | patient) + (1 | center)
##    Data: bp
## 
##      AIC      BIC   logLik deviance df.resid 
##   7492.4   7552.4  -3734.2   7468.4     1080 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9294 -0.5060 -0.0267  0.5195  3.4822 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  patient  (Intercept) 50.5730  7.1115        
##           visit        0.4893  0.6995   -0.55
##  center   (Intercept)  4.4599  2.1118        
##  Residual             33.0604  5.7498        
## Number of obs: 1092, groups:  patient, 288; center, 29
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            47.55432    8.89378 294.50121   5.347 1.80e-07 ***
## visit                  -0.57965    0.15608 276.84566  -3.714 0.000247 ***
## factor(treat)B         -1.59138    1.59385 276.19012  -0.998 0.318935    
## factor(treat)C         -3.26869    1.58166 273.63739  -2.067 0.039711 *  
## dbp1                    0.47959    0.08537 288.68482   5.618 4.55e-08 ***
## visit:factor(treat)B    0.05637    0.22121 272.92265   0.255 0.799035    
## visit:factor(treat)C    0.04252    0.21829 270.62461   0.195 0.845708    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) visit  fct()B fct()C dbp1   vs:()B
## visit       -0.098                                   
## factr(trt)B -0.123  0.554                            
## factr(trt)C -0.083  0.558  0.489                     
## dbp1        -0.991 -0.001  0.037 -0.005              
## vst:fctr()B  0.063 -0.706 -0.793 -0.394  0.007       
## vst:fctr()C  0.066 -0.715 -0.396 -0.794  0.005  0.505
  1. Compare the two models using the likelihood ratio test. Is the random intercept per center significant? Should you base your decision to include the intercept on this outcome?
anova(bp.2l.lme,bp.3l.lme)
## Data: bp
## Models:
## bp.2l.lme: dbp ~ visit + factor(treat) + visit * factor(treat) + dbp1 + 
## bp.2l.lme:     (visit | patient)
## bp.3l.lme: dbp ~ visit + factor(treat) + visit * factor(treat) + dbp1 + 
## bp.3l.lme:     (visit | patient) + (1 | center)
##           npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## bp.2l.lme   11 7501.6 7556.6 -3739.8   7479.6                         
## bp.3l.lme   12 7492.4 7552.4 -3734.2   7468.4 11.221  1  0.0008087 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding a random intercept per center significantly improves the model fit. It also reduces the variance of the patient random intercepts, indicating that some of the variation among patients is actually due to variation among centers. Since the design is 3-level, (at least) a random intercept per center should be included. The decision to include this random intercept per center should be based on the study design, and not on the results of testing the random intercept.

  1. Is the time*treatment interaction statistically significant? What does that mean for our answer to the research question?
bp.3l.lme2 <- lmer(dbp ~ visit + factor(treat)+ dbp1 + (visit|patient) + (1|center), data=bp, REML=F,
                   control = lmerControl(optimizer = "optimx", calc.derivs = FALSE, optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
## Loading required namespace: optimx
anova(bp.3l.lme, bp.3l.lme2)
## Data: bp
## Models:
## bp.3l.lme2: dbp ~ visit + factor(treat) + dbp1 + (visit | patient) + (1 | 
## bp.3l.lme2:     center)
## bp.3l.lme: dbp ~ visit + factor(treat) + visit * factor(treat) + dbp1 + 
## bp.3l.lme:     (visit | patient) + (1 | center)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## bp.3l.lme2   10 7488.5 7538.5 -3734.2   7468.5                     
## bp.3l.lme    12 7492.4 7552.4 -3734.2   7468.4 0.0707  2     0.9652

The treatment*time interaction can be dropped (LRT = 0.07 on 2 df, p almost 1)

Interpretation: starting from week 3, there is no difference in the patterns of DBP over time for the three treatment groups! That either means that there is no treatment effect in the study, or that the treatment effect was already present by week 3, the first post-randomization measurement we have here.

  1. Drop the interaction from the model, and test for a treatment effect.
# trying to drop treatment or dpb1 will result in significantly worse fits, this is the final model
bp.3l.lme3 <- lmer(dbp ~ visit + dbp1 + (visit|patient) + (1|center), data=bp, REML=F)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00282647 (tol = 0.002, component 1)
anova(bp.3l.lme2, bp.3l.lme3)
## Data: bp
## Models:
## bp.3l.lme3: dbp ~ visit + dbp1 + (visit | patient) + (1 | center)
## bp.3l.lme2: dbp ~ visit + factor(treat) + dbp1 + (visit | patient) + (1 | 
## bp.3l.lme2:     center)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## bp.3l.lme3    8 7494.3 7534.3 -3739.2   7478.3                        
## bp.3l.lme2   10 7488.5 7538.5 -3734.2   7468.5 9.8134  2   0.007397 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The model with fixed effects for treatment, time and baseline DBP is the final model. So there are significant treatment effect and time effects.

  1. Run the final model once more using REML estimation. See if you can get some post-hoc tests of treatments (hint: In R, the package lmerTest contains the useful functions ls_means() and difflsmeans().).
# run 1 more time, with REML
bp.3l.lmefin <- lmer(dbp ~ visit + factor(treat)+ dbp1 + (visit|patient) + (1|center), data=bp, REML=T)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00327907 (tol = 0.002, component 1)
summary(bp.3l.lmefin)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dbp ~ visit + factor(treat) + dbp1 + (visit | patient) + (1 |  
##     center)
##    Data: bp
## 
## REML criterion at convergence: 7470.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9044 -0.5108 -0.0248  0.5125  3.4560 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  patient  (Intercept) 51.4575  7.1734        
##           visit        0.4996  0.7068   -0.55
##  center   (Intercept)  4.7829  2.1870        
##  Residual             33.0481  5.7488        
## Number of obs: 1092, groups:  patient, 288; center, 29
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      47.48088    8.93239 287.22560   5.316 2.14e-07 ***
## visit            -0.54645    0.08976 268.98011  -6.088 3.93e-09 ***
## factor(treat)B   -1.27130    0.97662 265.65480  -1.302  0.19414    
## factor(treat)C   -3.02812    0.96637 263.31130  -3.134  0.00192 ** 
## dbp1              0.47860    0.08596 284.75942   5.568 5.97e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) visit  fct()B fct()C
## visit       -0.066                     
## factr(trt)B -0.122 -0.010              
## factr(trt)C -0.051 -0.019  0.487       
## dbp1        -0.994  0.009  0.071 -0.002
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00327907 (tol = 0.002, component 1)

Oops, difflsmeans(bp.3l.lmefin, test.effs='treat') does not work. It only works with factor variables, so first: make a new version of the treatment variable that is defined as a factor and use that in the model

bp$ftreat <- factor(bp$treat)
bp.3l.lmefin2 <- lmer(dbp ~ visit + ftreat + dbp1 + (visit|patient) + (1|center), data=bp, REML=T)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00327907 (tol = 0.002, component 1)
ls_means(bp.3l.lmefin2, test.effs="ftreat")
## Least Squares Means table:
## 
##          Estimate Std. Error   df t value    lower    upper  Pr(>|t|)    
## ftreatA  93.43826    0.83277 78.2  112.20 91.78042 95.09611 < 2.2e-16 ***
## ftreatB  92.16696    0.85353 81.4  107.98 90.46883 93.86510 < 2.2e-16 ***
## ftreatC  90.41014    0.83319 83.5  108.51 88.75310 92.06718 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Confidence level: 95%
##   Degrees of freedom method: Satterthwaite
difflsmeans(bp.3l.lmefin2, test.effs="ftreat")
## Least Squares Means table:
## 
##                     Estimate Std. Error    df t value    lower    upper
## ftreatA  - ftreatB   1.27130    0.97662 265.7  1.3017 -0.65161  3.19421
## ftreatA  - ftreatC   3.02812    0.96637 263.3  3.1335  1.12533  4.93092
## ftreatB  - ftreatC   1.75682    0.98372 260.8  1.7859 -0.18023  3.69387
##                     Pr(>|t|)   
## ftreatA  - ftreatB  0.194137   
## ftreatA  - ftreatC  0.001923 **
## ftreatB  - ftreatC  0.075278 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##   Confidence level: 95%
##   Degrees of freedom method: Satterthwaite

Bonferroni correction, p-values times the number of tests performed (e.g. 3) Confidence intervals, 1- α/(number of tests).

# a - b corrected 95%CI
qt(0.05/3, df=265.6, lower.tail = FALSE, log.p = FALSE)
## [1] 2.139175
1.3+(2.139175*0.977)
## [1] 3.389974
1.3-(2.139175*0.977)
## [1] -0.789974
# a - c corrected 95%CI
qt(0.05/3, df=263.3, lower.tail = FALSE, log.p = FALSE)
## [1] 2.139273
3.0-(2.139273*0.966)
## [1] 0.9334623
3.0+(2.139273*0.966)
## [1] 5.066538
# b - b corrected 95%CI
qt(0.05/3, df=260.8, lower.tail = FALSE, log.p = FALSE)
## [1] 2.139381
1.8-(2.139381*0.984)
## [1] -0.3051509
1.8+(2.138381*0.984)
## [1] 3.904167
  1. Interpret the results in the context of the research question.

When we control for treatment and baseline DBP, patients’ DBP decreases by a little over half a mmHg per week on average; this decrease is estimated to be the same for all three treatment groups (because the treat*visit interaction was not significant and we therefore left it out of the model). There is a significant treatment effect: after adjustment for multiple comparison, treatment C reduces blood pressure compared to treatment A by 3 points.

Exercise 5

Centering the explanatory variables is a major topic in mixed models. (See for example the evaluation review of Omar Paccagnella “Centering or Not Centering in Multilevel Models? The Role of the Group Mean and the Assessment of Group Effects”.) In this exercise we will demonstrate the consequences of centering. Consider again the longitudinal data example of the paper by John Fox (see exercise 2 of assignment part 2), available as Blackmore in the car package.

library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:lme4':
## 
##     lmList
data(Blackmore)
Blackmore$log.exercise <- log(Blackmore$exercise + 5/60, 2)
  1. Fit a model with random intercept, random age effect, fixed group and a fixed interaction between age and group. Don’t forget the transformation of the outcome! Make a table of the most important parameter estimates.
bm.lme.a <- lme(log.exercise ~ age*group, random=~age|subject, data=Blackmore)
summary(bm.lme.a)
## Linear mixed-effects model fit by REML
##  Data: Blackmore 
##        AIC      BIC    logLik
##   3630.136 3668.912 -1807.068
## 
## Random effects:
##  Formula: ~age | subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 2.2115383 (Intr)
## age         0.1647956 -0.779
## Residual    1.2440950       
## 
## Fixed effects: log.exercise ~ age * group 
##                       Value Std.Error  DF   t-value p-value
## (Intercept)      -0.7881948 0.3754620 712 -2.099266  0.0361
## age               0.0640222 0.0313605 712  2.041491  0.0416
## grouppatient     -2.2728627 0.4768124 229 -4.766786  0.0000
## age:grouppatient  0.2398586 0.0394074 712  6.086645  0.0000
##  Correlation: 
##                  (Intr) age    grpptn
## age              -0.906              
## grouppatient     -0.787  0.713       
## age:grouppatient  0.721 -0.796 -0.903
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7348640 -0.4245452  0.1227672  0.5280055  2.6361915 
## 
## Number of Observations: 945
## Number of Groups: 231
compareCoefs(bm.lme.a)
## Calls:
## lme.formula(fixed = log.exercise ~ age * group, data = Blackmore, random = 
##   ~age | subject)
## 
##                  Model 1
## (Intercept)       -0.788
## SE                 0.375
##                         
## age               0.0640
## SE                0.0314
##                         
## grouppatient      -2.273
## SE                 0.477
##                         
## age:grouppatient  0.2399
## SE                0.0394
##
  1. Repeat the analysis but first transform the age to age-8. Extend the table of (a) with the new estimates.
bm.lme.b <- lme(log.exercise ~ I(age-8)*group, random = ~I(age-8)|subject, data=Blackmore)
summary(bm.lme.b)
## Linear mixed-effects model fit by REML
##  Data: Blackmore 
##        AIC      BIC    logLik
##   3630.136 3668.912 -1807.068
## 
## Random effects:
##  Formula: ~I(age - 8) | subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 1.4435580 (Intr)
## I(age - 8)  0.1647954 -0.281
## Residual    1.2440951       
## 
## Fixed effects: log.exercise ~ I(age - 8) * group 
##                              Value  Std.Error  DF   t-value p-value
## (Intercept)             -0.2760170 0.18236870 712 -1.513511  0.1306
## I(age - 8)               0.0640222 0.03136052 712  2.041491  0.0416
## grouppatient            -0.3539943 0.23529124 229 -1.504494  0.1338
## I(age - 8):grouppatient  0.2398585 0.03940734 712  6.086646  0.0000
##  Correlation: 
##                         (Intr) I(g-8) grpptn
## I(age - 8)              -0.489              
## grouppatient            -0.775  0.379       
## I(age - 8):grouppatient  0.389 -0.796 -0.489
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7348639 -0.4245457  0.1227672  0.5280054  2.6361916 
## 
## Number of Observations: 945
## Number of Groups: 231
compareCoefs(bm.lme.a,bm.lme.b)
## Calls:
## 1: lme.formula(fixed = log.exercise ~ age * group, data = Blackmore, random 
##   = ~age | subject)
## 2: lme.formula(fixed = log.exercise ~ I(age - 8) * group, data = Blackmore, 
##   random = ~I(age - 8) | subject)
## 
##                         Model 1 Model 2
## (Intercept)              -0.788  -0.276
## SE                        0.375   0.182
##                                        
## age                      0.0640        
## SE                       0.0314        
##                                        
## grouppatient             -2.273  -0.354
## SE                        0.477   0.235
##                                        
## age:grouppatient         0.2399        
## SE                       0.0394        
##                                        
## I(age - 8)                       0.0640
## SE                               0.0314
##                                        
## I(age - 8):grouppatient          0.2399
## SE                               0.0394
##
  1. Repeat the analysis again but now transform the age to the deviation from the mean age (i.e. centering to the mean). Extend again the table with the new estimates (recall the compareCoefs() function).
bm.lme.c <- lme(log.exercise ~ I(age-mean(age))*group, random = ~I(age-mean(age))|subject, data=Blackmore)
summary(bm.lme.c)
## Linear mixed-effects model fit by REML
##  Data: Blackmore 
##        AIC      BIC    logLik
##   3630.136 3668.912 -1807.068
## 
## Random effects:
##  Formula: ~I(age - mean(age)) | subject
##  Structure: General positive-definite, Log-Cholesky parametrization
##                    StdDev    Corr  
## (Intercept)        1.3948727 (Intr)
## I(age - mean(age)) 0.1647955 0.116 
## Residual           1.2440951       
## 
## Fixed effects: log.exercise ~ I(age - mean(age)) * group 
##                                      Value  Std.Error  DF   t-value p-value
## (Intercept)                     -0.0556742 0.16014718 712 -0.347644  0.7282
## I(age - mean(age))               0.0640222 0.03136052 712  2.041491  0.0416
## grouppatient                     0.4715176 0.20621237 229  2.286563  0.0231
## I(age - mean(age)):grouppatient  0.2398585 0.03940735 712  6.086645  0.0000
##  Correlation: 
##                                 (Intr) I(-mn()) grpptn
## I(age - mean(age))               0.117                
## grouppatient                    -0.777 -0.091         
## I(age - mean(age)):grouppatient -0.093 -0.796    0.099
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7348639 -0.4245456  0.1227672  0.5280054  2.6361916 
## 
## Number of Observations: 945
## Number of Groups: 231
compareCoefs(bm.lme.a,bm.lme.b,bm.lme.c)
## Calls:
## 1: lme.formula(fixed = log.exercise ~ age * group, data = Blackmore, random 
##   = ~age | subject)
## 2: lme.formula(fixed = log.exercise ~ I(age - 8) * group, data = Blackmore, 
##   random = ~I(age - 8) | subject)
## 3: lme.formula(fixed = log.exercise ~ I(age - mean(age)) * group, data = 
##   Blackmore, random = ~I(age - mean(age)) | subject)
## 
##                                 Model 1 Model 2 Model 3
## (Intercept)                     -0.7882 -0.2760 -0.0557
## SE                               0.3755  0.1824  0.1601
##                                                        
## age                              0.0640                
## SE                               0.0314                
##                                                        
## grouppatient                     -2.273  -0.354   0.472
## SE                                0.477   0.235   0.206
##                                                        
## age:grouppatient                 0.2399                
## SE                               0.0394                
##                                                        
## I(age - 8)                               0.0640        
## SE                                       0.0314        
##                                                        
## I(age - 8):grouppatient                  0.2399        
## SE                                       0.0394        
##                                                        
## I(age - mean(age))                               0.0640
## SE                                               0.0314
##                                                        
## I(age - mean(age)):grouppatient                  0.2399
## SE                                               0.0394
##
  1. Compare the results of these three models and explain the differences and similarities (Hint. Use something like figure 9 of the Fox & Weisberg paper).
# Plot predictions for model b
plot(Blackmore$age, predict(bm.lme.b, level=0))

# Check residuals
plot(bm.lme.a)

### plots to compare models
xfrom0= seq(0,10,length.out=100)
xfrom8= seq(0,18, length.out=100)
xcenter= Blackmore$age-mean(Blackmore$age)
par(mfrow=c(1,3))
plot(xfrom0, -0.06386+(xfrom0*0.014), type="l", ylim=c(-0.8,0.8), xlab="x starting at age 8, 0=8", ylab="predicted log.exercise")
points(xfrom0, -0.06386+(xfrom0*0.014)+(xfrom0*0.0555) -0.081 , type="l", col="red")
plot(xfrom8, -0.18+(xfrom8*0.014), type="l", ylim=c(-0.8,0.8), xlab="x starting at age 0", xlim=c(0,18), ylab="predicted log.exercise")
points(xfrom8, -0.18+(xfrom8*0.014)+(xfrom8*0.0555) -0.52589 , type="l", col="red")
plot(xcenter, -0.01288+(xcenter*0.00726), type="l", ylim=c(-0.8,0.8), xlab="age centered", xlim=range(xcenter), ylab="predicted log.exercise")
points(xcenter, -0.01288+(xcenter*0.00726) + (xcenter*0.055498) +0.1091 , type="l", col="red")

par(mfrow=c(1,1))
  1. Which model would you prefer? Explain your choice.

As we have already seen with the Reisby dataset, the choice of a “0” value for time (age, in this case) changes the value and interpretation of the fixed intercept, along with the variance of the random intercepts and the correlation of the random intercepts and slopes per girl. Since the girls in the study are at least 8 at the beginning of the study, having age start at 0 is a bit of an odd choice; the intercept (and variation around it) estimates the mean energy expenditure of the girls as newborns and the estimated variation of the random intercepts is also being based on a time point well outside the range of the data. Choosing age-8 means the null point is at age 8, so the fixed intercept is the estimated mean (log) energy expenditure at age 8 and the random effect indicates the variation around that. Centering age gives us a null point of 11.4 years, and the intercept will then be the estimated mean (log) energy expenditure at age 11.4.

Note that the estimated main effect of age and the interaction of age and group are unaffected by the change in parameterization of age, but the group effect is quite different for the 3 models. The estimate for group is the mean difference in (log) energy expenditure for the patients vs. controls at ages 0, 8 or 11.4. Since the group*age interaction is significant, the difference between the two groups depends on age.

Exercise 6 (CHALLENGE EXERCISE)

Kroesbergen et al. (Eur Respir J 1999) investigated the flow-dependency of exhaled nitric oxide (NO) in healthy children (n=20) and children with asthma (n=19). The concentration of NO in exhaled are was measured four times, at four different target values of the flow of exhalation (2, 5, 10, and 20% of their vital capacity per second). The actual flows differed from the target flows.

The following questions are addressed: is there is an association between NO (pbb) and FLOW (L/sec) for healthy and asthmatic children? Is this association different for the two groups? The variables NO and FLOW have been been log-transformed.

The dataset no.Rdata contains the following variables:

  • NUM: child’s identification number
  • DIAGNOSE: 0 = “healthy” and 1 = “asthma”
  • FLOW: 10log(flow)-10log(2)
  • NO: 10log(NO concentration)
# read in data
load(file = "no.Rdata")
# Make a groupedData object for the "trellis graphs" (spaghetti plot) per group
nogr <- nlme::groupedData(no ~ flow | num, data=nodat, outer = ~diagnose)
# spaghetti plot per group
plot(nogr, outer=~diagnose, key=F, aspect=1)

  1. Describe what you see in the graphs. What would be your first intuition with respect to the research questions?

There seems to be a negative association between NO and flow. Clearly random intercepts are necessary for both groups, possibly also random slopes. The slopes seem to be steeper in the asthmatic group. There is also more variation between children in the asthmatic group than in the healthy group.

  1. Fit a model with random intercept random slope regression for the healthy group. What is the mean slope and intercept? What are the estimated variances of the random components? Is a random slope model necessary, or would a random intercept be sufficient?
# make subsets of healthy and asthmatic children
healthy <- subset(nodat,nodat$diagnose==0)
asthma <- subset(nodat,nodat$diagnose==1)
# mixed model: NO as function of flow, random int+slope, healthy group
no.lme.h <- lme(no~flow , random=~flow| num, method="ML", data=healthy)
summary(no.lme.h)
## Linear mixed-effects model fit by maximum likelihood
##  Data: healthy 
##         AIC       BIC   logLik
##   -162.0047 -147.7125 87.00234
## 
## Random effects:
##  Formula: ~flow | num
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 0.1372792 (Intr)
## flow        0.1018528 -0.839
## Residual    0.0496389       
## 
## Fixed effects: no ~ flow 
##                  Value  Std.Error DF   t-value p-value
## (Intercept)  0.8108212 0.03162558 59  25.63815       0
## flow        -0.3452360 0.02754635 59 -12.53291       0
##  Correlation: 
##      (Intr)
## flow -0.707
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.15285664 -0.58959632 -0.03373757  0.48607120  2.46376179 
## 
## Number of Observations: 80
## Number of Groups: 20

Mean slope for healthy children is -0.345 (se 0.0275), mean intercept 0.811. Standard deviation for random intercept is 0.137 (variance = 0.0188), and for the random slope 0.102 (var = 0.0104).

# check random slope for healthy group
no.lme.h2 <- update(no.lme.h, random= ~1|num)
anova(no.lme.h, no.lme.h2)
##           Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## no.lme.h      1  6 -162.0047 -147.7125 87.00234                        
## no.lme.h2     2  4 -143.5307 -134.0026 75.76534 1 vs 2 22.47401  <.0001

The random slope is necessary.

  1. Answer the same questions for the group of children with asthma.
# mixed model: NO as function of flow, random int+slope, asthmatic group
no.lme.a <- lme(no~flow , random=~flow| num, method="ML", data=asthma)
summary(no.lme.a)
## Linear mixed-effects model fit by maximum likelihood
##  Data: asthma 
##        AIC      BIC  logLik
##   -64.1888 -50.2044 38.0944
## 
## Random effects:
##  Formula: ~flow | num
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept) 0.38293951 (Intr)
## flow        0.18164821 -0.393
## Residual    0.06682339       
## 
## Fixed effects: no ~ flow 
##                  Value  Std.Error DF    t-value p-value
## (Intercept)  0.8517394 0.08940237 56   9.527034       0
## flow        -0.5848256 0.04793471 56 -12.200462       0
##  Correlation: 
##      (Intr)
## flow -0.346
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.61299643 -0.53347738  0.04039453  0.38132723  2.06154236 
## 
## Number of Observations: 76
## Number of Groups: 19
# check random slope for asthma group
no.lme.a2 <- update(no.lme.a, random= ~1|num)
summary(no.lme.a2)
## Linear mixed-effects model fit by maximum likelihood
##  Data: asthma 
##         AIC       BIC   logLik
##   -46.52904 -37.20611 27.26452
## 
## Random effects:
##  Formula: ~1 | num
##         (Intercept)  Residual
## StdDev:   0.3762181 0.1024387
## 
## Fixed effects: no ~ flow 
##                  Value  Std.Error DF    t-value p-value
## (Intercept)  0.8513291 0.08827587 56   9.643962       0
## flow        -0.5776986 0.03358039 56 -17.203453       0
##  Correlation: 
##      (Intr)
## flow -0.001
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.731791522 -0.497431958  0.008702399  0.509985217  2.035226345 
## 
## Number of Observations: 76
## Number of Groups: 19
anova(no.lme.a, no.lme.a2)
##           Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## no.lme.a      1  6 -64.18880 -50.20440 38.09440                        
## no.lme.a2     2  4 -46.52904 -37.20611 27.26452 1 vs 2 21.65976  <.0001

Mean slope for asthmatic children is -0.585 (se 0.0479), mean intercept 0.852. Standard deviation for the random intercept is 0.383 (variance = 0.147) and for random slope 0.182 (var=0.0330). The random slope is also necessary here.

  1. Compare the results of the two groups. What would you conclude on the basis of these results?

The mean slope appears to be different (-0.345 for healthy vs -0.585 for asthmatic children) for the two groups with asthmatics exhaling less NO at the higher flow rates. The variation in the random intercepts and slopes is larger for the asthmatic children.

  1. Fit a model for the two groups combined, using an interaction for diagnosis*flow and allowing a different covariance structure for the two groups (note: In R, use the following command:)
no.lme.3 <- lme(no~flow + factor(diagnose)+ flow:factor(diagnose), random=~flow + diagnose + flow*diagnose| num,  weights=varIdent(form=~1|diagnose), method="ML", data=nodat)

The random command contains a term for group, flow and the group*flow interaction and thus allows for a separate variance for the random intercept and slope for the two diagnosis groups. The weights statement further allows differing residual variances for the two groups. Check that the log-likelihood of this model is the sum of the log-likelihoods of the models in b and c.

summary(no.lme.3)
## Linear mixed-effects model fit by maximum likelihood
##  Data: nodat 
##         AIC       BIC   logLik
##   -218.1935 -169.3958 125.0967
## 
## Random effects:
##  Formula: ~flow + diagnose + flow * diagnose | num
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev     Corr                
## (Intercept)   0.13727910 (Intr) flow   diagns
## flow          0.10185272 -0.839              
## diagnose      0.26189009  0.823 -0.680       
## flow:diagnose 0.17568175 -0.014 -0.230  0.062
## Residual      0.04963893                     
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | diagnose 
##  Parameter estimates:
##        0        1 
## 1.000000 1.346191 
## Fixed effects: no ~ flow + factor(diagnose) + flow:factor(diagnose) 
##                             Value  Std.Error  DF    t-value p-value
## (Intercept)             0.8108212 0.03163596 115  25.629735  0.0000
## flow                   -0.3452360 0.02755540 115 -12.528798  0.0000
## factor(diagnose)1       0.0409182 0.09480556  37   0.431601  0.6685
## flow:factor(diagnose)1 -0.2395895 0.05527596 115  -4.334425  0.0000
##  Correlation: 
##                        (Intr) flow   fct()1
## flow                   -0.707              
## factor(diagnose)1      -0.334  0.236       
## flow:factor(diagnose)1  0.353 -0.499 -0.400
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.612996113 -0.565060550  0.009625799  0.435834542  2.463760385 
## 
## Number of Observations: 156
## Number of Groups: 39

Note that there are now 4 different standard deviations for the random effects: one each for intercept and slope for the healthy group (0.137 and 0.102, respectively) and one each for intercept and slope for the asthma group (unfortunately, these cannot be read easily from the output; with variances we cannot simply add or subtract main effects and interactions as we can with fixed effects). Note that the standard deviations for the healthy group are the same as in the stratified analysis in (b). Note also that the log-likelihood of the full model in (b) was 87.002 and that the log-likelihood of the full model in (c) was 38.094. If we add these together, we get 125.09, which is the log-likelihood of the current model. The current model is thus another way of doing a completely stratified model, as in (b) and (c).

  1. Can the variance structure be simplified by assuming equal residual variances? (Hint: try removing the weights statement.)
# reduce variance structure: take out differing residual variances
no.lme.4 <- lme(no~flow + factor(diagnose)+ flow:factor(diagnose), 
                random=~flow + diagnose + flow*diagnose| num, method="ML", data=nodat)
summary(no.lme.4)
## Linear mixed-effects model fit by maximum likelihood
##  Data: nodat 
##         AIC       BIC   logLik
##   -216.8824 -171.1346 123.4412
## 
## Random effects:
##  Formula: ~flow + diagnose + flow * diagnose | num
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev     Corr                
## (Intercept)   0.13659185 (Intr) flow   diagns
## flow          0.08989251 -0.938              
## diagnose      0.26884370  0.765 -0.749       
## flow:diagnose 0.17504544  0.116 -0.128 -0.005
## Residual      0.05896974                     
## 
## Fixed effects: no ~ flow + factor(diagnose) + flow:factor(diagnose) 
##                             Value  Std.Error  DF    t-value p-value
## (Intercept)             0.8113303 0.03168699 115  25.604525  0.0000
## flow                   -0.3436318 0.02705552 115 -12.700987  0.0000
## factor(diagnose)1       0.0404842 0.09488772  37   0.426653  0.6721
## flow:factor(diagnose)1 -0.2417015 0.05485093 115  -4.406517  0.0000
##  Correlation: 
##                        (Intr) flow   fct()1
## flow                   -0.713              
## factor(diagnose)1      -0.334  0.238       
## flow:factor(diagnose)1  0.352 -0.493 -0.400
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.91489213 -0.53349386 -0.01647227  0.47910133  2.36872501 
## 
## Number of Observations: 156
## Number of Groups: 39
# is it necessary to allow for differing residual variances?
anova(no.lme.3,no.lme.4)
##          Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## no.lme.3     1 16 -218.1935 -169.3958 125.0967                        
## no.lme.4     2 15 -216.8824 -171.1346 123.4412 1 vs 2 3.311073  0.0688

p = 0.0688 / 2 = 0.0344. It appears the allowing for differing residual variances is necessary.

  1. Can the covariance structure be further reduced? Try to simplify the random command by removing the diagnose*flow interaction, and then by removing diagnose.
# it would appear that the residual variances differ. Leave weights statement in.
# is it necessary to have different random slope variances for the two groups?
no.lme.5 <- lme(no~flow + factor(diagnose)+ flow:factor(diagnose), 
                random=~flow + diagnose| num, weights=varIdent(form=~1|diagnose),
                method="ML", data=nodat)
anova(no.lme.3,no.lme.5)
##          Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## no.lme.3     1 16 -218.1935 -169.3958 125.0967                        
## no.lme.5     2 12 -222.7014 -186.1031 123.3507 1 vs 2 3.492107  0.4791

p = 0.4791 / 2 = 0.2396, separate variances for random slope not necessary.

# try to reduce variance structure further: necessary to have different
# random intercept variances for the two groups?
no.lme.6 <- lme(no~flow + factor(diagnose)+ flow:factor(diagnose), 
                random=~flow | num, weights=varIdent(form=~1|diagnose),
                method="ML", data=nodat)
anova(no.lme.5,no.lme.6)
##          Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## no.lme.5     1 12 -222.7014 -186.1031 123.3507                        
## no.lme.6     2  9 -203.5273 -176.0786 110.7636 1 vs 2 25.17408  <.0001

p < 0.0001 / 2 = 0.00005, separate variances for random intercept necessary.

  1. Is the fixed diagnose*flow interaction necessary? In other words: do the asthmatic children indeed have a significantly more negative slope than the healthy children?
# is the fixed flow*group interaction necessary? model without interaction
no.lme.7 <- lme(no~flow + factor(diagnose), random=~flow + diagnose| num,
                weights=varIdent(form=~1|diagnose), method="ML", data=nodat)
anova(no.lme.5,no.lme.7)
##          Model df       AIC       BIC   logLik   Test L.Ratio p-value
## no.lme.5     1 12 -222.7014 -186.1031 123.3507                       
## no.lme.7     2 11 -208.2011 -174.6527 115.1005 1 vs 2 16.5003  <.0001

p < 0.0001, the interaction is significant. The asthmatic children have a significantly more negative slope than the healthy children.

  1. Answer the following questions for the model you end up with. What is the slope for the healthy group, what for the asthma group? What is the standard deviation of the random intercept the healthy children, and is the variance for asthmatic children larger or smaller? What is the residual standard deviation for the two groups?
# Final model (#5), with REML
no.lme.final <- lme(no~flow + factor(diagnose)+ flow:factor(diagnose), 
                    random=~flow + diagnose| num, weights=varIdent(form=~1|diagnose),
                    method="REML", data=nodat)
summary(no.lme.final)
## Linear mixed-effects model fit by REML
##  Data: nodat 
##         AIC      BIC   logLik
##   -204.3626 -168.076 114.1813
## 
## Random effects:
##  Formula: ~flow + diagnose | num
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr         
## (Intercept) 0.16080858 (Intr) flow  
## flow        0.14515225 -0.852       
## diagnose    0.29348041  0.384  0.001
## Residual    0.04850413              
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | diagnose 
##  Parameter estimates:
##        0        1 
## 1.000000 1.420661 
## Fixed effects: no ~ flow + factor(diagnose) + flow:factor(diagnose) 
##                             Value  Std.Error  DF   t-value p-value
## (Intercept)             0.8103540 0.03639792 115 22.263748  0.0000
## flow                   -0.3464639 0.03557681 115 -9.738475  0.0000
## factor(diagnose)1       0.0413823 0.09587913  37  0.431609  0.6685
## flow:factor(diagnose)1 -0.2374768 0.05385543 115 -4.409524  0.0000
##  Correlation: 
##                        (Intr) flow   fct()1
## flow                   -0.778              
## factor(diagnose)1      -0.380  0.295       
## flow:factor(diagnose)1  0.514 -0.661 -0.398
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.58783073 -0.53427169  0.04976964  0.47714277  2.26198962 
## 
## Number of Observations: 156
## Number of Groups: 39

Mean intercept and slope for healthy children: 0.810 and -0.346. Mean intercept and slope for asthmatic children: (0.810+0.041=) 0.852 and (-0.346-0.237=) -0.584.

Standard deviation random intercept healthy children: 0.161 (var = 0.0259). Standard deviation random intercept asthmatic children cannot easily be calculated, but is larger than that of healthy children.

Residual standard deviation for healthy children is 0.048; for asthmatic children 1.42*0.0485 = 0.0689

Session Info

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.3 (2020-10-10)
##  os       macOS Catalina 10.15.7      
##  system   x86_64, darwin17.0          
##  ui       X11                         
##  language (EN)                        
##  collate  nl_NL.UTF-8                 
##  ctype    nl_NL.UTF-8                 
##  tz       Europe/Amsterdam            
##  date     2020-11-30                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version    date       lib source        
##  abind         1.4-5      2016-07-21 [1] CRAN (R 4.0.2)
##  assertthat    0.2.1      2019-03-21 [1] CRAN (R 4.0.2)
##  boot          1.3-25     2020-04-26 [2] CRAN (R 4.0.3)
##  callr         3.5.1      2020-10-13 [1] CRAN (R 4.0.2)
##  car         * 3.0-10     2020-09-29 [1] CRAN (R 4.0.2)
##  carData     * 3.0-4      2020-05-22 [1] CRAN (R 4.0.2)
##  cellranger    1.1.0      2016-07-27 [1] CRAN (R 4.0.2)
##  cli           2.2.0      2020-11-20 [1] CRAN (R 4.0.2)
##  colorspace    2.0-0      2020-11-11 [1] CRAN (R 4.0.2)
##  crayon        1.3.4      2017-09-16 [1] CRAN (R 4.0.2)
##  curl          4.3        2019-12-02 [1] CRAN (R 4.0.1)
##  data.table    1.13.2     2020-10-19 [1] CRAN (R 4.0.2)
##  desc          1.2.0      2018-05-01 [1] CRAN (R 4.0.2)
##  devtools      2.3.2      2020-09-18 [1] CRAN (R 4.0.2)
##  digest        0.6.27     2020-10-24 [1] CRAN (R 4.0.2)
##  dplyr         1.0.2      2020-08-18 [1] CRAN (R 4.0.2)
##  ellipsis      0.3.1      2020-05-15 [1] CRAN (R 4.0.2)
##  evaluate      0.14       2019-05-28 [1] CRAN (R 4.0.1)
##  fansi         0.4.1      2020-01-08 [1] CRAN (R 4.0.2)
##  forcats       0.5.0      2020-03-01 [1] CRAN (R 4.0.2)
##  foreign       0.8-80     2020-05-24 [2] CRAN (R 4.0.3)
##  fs            1.5.0      2020-07-31 [1] CRAN (R 4.0.2)
##  generics      0.1.0      2020-10-31 [1] CRAN (R 4.0.2)
##  ggplot2     * 3.3.2      2020-06-19 [1] CRAN (R 4.0.2)
##  glue          1.4.2      2020-08-27 [1] CRAN (R 4.0.2)
##  gtable        0.3.0      2019-03-25 [1] CRAN (R 4.0.2)
##  haven         2.3.1      2020-06-01 [1] CRAN (R 4.0.2)
##  hms           0.5.3      2020-01-08 [1] CRAN (R 4.0.2)
##  htmltools     0.5.0      2020-06-16 [1] CRAN (R 4.0.2)
##  knitr         1.30       2020-09-22 [1] CRAN (R 4.0.2)
##  lattice       0.20-41    2020-04-02 [2] CRAN (R 4.0.3)
##  lifecycle     0.2.0      2020-03-06 [1] CRAN (R 4.0.2)
##  lme4        * 1.1-25     2020-10-23 [1] CRAN (R 4.0.2)
##  lmerTest    * 3.1-3      2020-10-23 [1] CRAN (R 4.0.2)
##  magrittr      2.0.1      2020-11-17 [1] CRAN (R 4.0.3)
##  MASS          7.3-53     2020-09-09 [2] CRAN (R 4.0.3)
##  Matrix      * 1.2-18     2019-11-27 [2] CRAN (R 4.0.3)
##  memoise       1.1.0      2017-04-21 [1] CRAN (R 4.0.2)
##  minqa         1.2.4      2014-10-09 [1] CRAN (R 4.0.2)
##  munsell       0.5.0      2018-06-12 [1] CRAN (R 4.0.2)
##  nlme        * 3.1-149    2020-08-23 [2] CRAN (R 4.0.3)
##  nloptr        1.2.2.2    2020-07-02 [1] CRAN (R 4.0.2)
##  numDeriv      2016.8-1.1 2019-06-06 [1] CRAN (R 4.0.2)
##  openxlsx      4.2.3      2020-10-27 [1] CRAN (R 4.0.2)
##  optimx        2020-4.2   2020-04-08 [1] CRAN (R 4.0.2)
##  pillar        1.4.7      2020-11-20 [1] CRAN (R 4.0.2)
##  pkgbuild      1.1.0      2020-07-13 [1] CRAN (R 4.0.2)
##  pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 4.0.2)
##  pkgload       1.1.0      2020-05-29 [1] CRAN (R 4.0.2)
##  prettyunits   1.1.1      2020-01-24 [1] CRAN (R 4.0.2)
##  processx      3.4.4      2020-09-03 [1] CRAN (R 4.0.2)
##  ps            1.4.0      2020-10-07 [1] CRAN (R 4.0.2)
##  purrr         0.3.4      2020-04-17 [1] CRAN (R 4.0.2)
##  R6            2.5.0      2020-10-28 [1] CRAN (R 4.0.2)
##  Rcpp          1.0.5      2020-07-06 [1] CRAN (R 4.0.2)
##  readxl        1.3.1      2019-03-13 [1] CRAN (R 4.0.2)
##  remotes       2.2.0      2020-07-21 [1] CRAN (R 4.0.2)
##  rio           0.5.16     2018-11-26 [1] CRAN (R 4.0.2)
##  rlang         0.4.8      2020-10-08 [1] CRAN (R 4.0.2)
##  rmarkdown     2.5        2020-10-21 [1] CRAN (R 4.0.3)
##  rprojroot     2.0.2      2020-11-15 [1] CRAN (R 4.0.2)
##  scales        1.1.1      2020-05-11 [1] CRAN (R 4.0.2)
##  sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 4.0.2)
##  statmod       1.4.35     2020-10-19 [1] CRAN (R 4.0.2)
##  stringi       1.5.3      2020-09-09 [1] CRAN (R 4.0.2)
##  stringr       1.4.0      2019-02-10 [1] CRAN (R 4.0.2)
##  testthat      3.0.0      2020-10-31 [1] CRAN (R 4.0.2)
##  tibble        3.0.4      2020-10-12 [1] CRAN (R 4.0.2)
##  tidyselect    1.1.0      2020-05-11 [1] CRAN (R 4.0.2)
##  usethis       1.6.3      2020-09-17 [1] CRAN (R 4.0.2)
##  vctrs         0.3.5      2020-11-17 [1] CRAN (R 4.0.2)
##  withr         2.3.0      2020-09-22 [1] CRAN (R 4.0.2)
##  xfun          0.19       2020-10-30 [1] CRAN (R 4.0.2)
##  yaml          2.2.1      2020-02-01 [1] CRAN (R 4.0.2)
##  zip           2.1.1      2020-08-27 [1] CRAN (R 4.0.2)
## 
## [1] /Users/dveen5/Library/R/4.0/library
## [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library